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Abstract 



The fully nonlinear and weakly dispersive Green-Naghdi model for shallow wa- 
ter waves of large amplitude is studied. The original model is first recast under 
r^ , a new formulation more suitable for numerical resolution. An hybrid finite vol- 

ume and finite difference splitting approach is then proposed. The hyperbolic 
part of the equations is handled with a high-order finite volume scheme allowing 
for breaking waves and dry areas. The dispersive part is treated with a classical 
finite difference approach. Extensive numerical validations are then performed 
in one horizontal dimension, relying both on analytical solutions and experi- 
mental data. The results show that our approach gives a good account of all 
the processes of wave transformation in coastal areas: shoaling, wave breaking 



iy-> , and run-up. 

Keywords: Green-Naghdi model, nonlinear shallow water, splitting method, 
^C^ \ finite volume, high order relaxation scheme, run-up. 
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1. Introduction 

^^ , In an incompressible, homogeneous, inviscid fluid, the propagation of surface 

$H ' waves is governed by the Euler equations with nonlinear boundary conditions 

at the surface and at the bottom. In its full generality, this problem is very 
complicated to solve, both mathematically and numerically. This is the reason 
why more simple models have been derived to describe the behavior of the 
solution in some physical specific regimes. A recent review of the different 
models that can be derived can be found in [2J|. 
Of particular interest in coastal oceanography is the shallow-water regime, which 
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corresponds to the configuration where the wave length A of the flow is large 
compared to the typical depth ho'. 

(Shallow water regime) /^ := — 2- <c 1. 

When the typical amplitude a of the wave is small, in the sense that 

(Small amplitude regime) e := — = 0{i-i), 

ho 

it is known that an approximation of order 0{y?) of the free surface Euler 
equations is furnished by the Boussinesq systems, such as the one derived by 



Peregrine in [32| for uneven bottoms. This model couples the surface elevation 
C, to the vertically averaged horizontal component of the velocity V, and can be 
written in non-dimensionalized form as 

dtV + e{V ■V)V + VC,= ^iD + 0{^l'^), ^' 

where h is the water depth and T) accounts for the nonhydrostatic and dispersive 
effects, and is a function of C, V^ and their derivatives. For instance, in the 



Boussinesq model derived in [S^ , one has 

V = ^V[V • {hdtV)] - ^V^dtV. (2) 

Unfortunately, the small amplitude assumption e ~ 0(/i) is too restrictive for 
many applications in coastal oceanography, where large amplitude waves have 
to be considered, 

(Large amplitude regime) e := -— = 0(1). 

ho 

If one wants to keep the same 0{iJ?) precision of ([1]) in a large amplitude regime, 
then the expression for T) is much more complicated than in ^. For instance, 
in II? and for flat bottoms, one has 

v - ^d.Ah'^iV^t + vv,., - (K)')]. 

In this regime, the corresponding equations ([l} have been derived flrst by Serre 



and then Su and Gardner [35| , Seabra-Santos et al. 3J] and Green and Naghdi 
[19| (other relevant references are [ij, |40|, |31|); consequently, these equations 
carry several names: Serre, Green-Naghdi, or fully nonlinear Boussinesq equa- 
tions. We will call them Green-Naghdi equations throughout this paper. Here 
again, we refer to [2J] for more details; note also that a rigorous mathematical 
justification of these models has been given in [l|. 

The Green-Naghdi equations ([T]) provide a correct description of the waves 
up to the breaking point; from this point however, they become useless (at least 



without consequent modifications). A first approach to model wave breaking 
is to add an ad hoc viscous term to the momentum equation, whose role is to 
account for the energy dissipation that occurs during wave breaking. This ap- 
proach has been used for instance by Zelt [4j] or Kennedy [22[ and Chen |13j . 
Recently, Cienfuegos et al. [16| proposed a new ID wave-breaking parametriza- 
tion including viscous-like effects on both the mass and the momentum equa- 
tions. This approach is able to reproduce wave height decay and intraphase 
nonlinear properties within the entire surf zone. However, the extension of this 
ad hoc parametrization to 2D wave cases remains a very difficult task. Another 
approach to handle wave breaking is to use the classical nonlinear shallow water 
equations, defined with I? = in ([!]) and denoted by NSWE in the follow- 
ing. These equations being hyperbolic, they develop shocks; after the breaking 
point, the waves are then described by the weak solutions of this hyperbolic 
system. This approach, used in [23] and [g is satisfactory in the sense that it 
gives a natural and correct description of the dissipation of energy during wave 
breaking. Its drawback, however, is that it is inappropriate in the shoaling zone 
since this models neglects the nonhydrostatic and dispersive effects. The moti- 
vation of this paper is to develop a model and a numerical scheme that describes 
correctly both phenomena. More precisely, we want to 

tt 1 Provide a good description of the dispersive effects (in the shoaling zone 

in particular); 
U 2 Take into account wave breaking in a simple way. 

Another theoretical and numerical difficulty in coastal oceanography is the 
description of the shoreline, i.e. the zone where the water depth vanishes, as 
the size of the computational domain becomes part of the solution. Taking 
into account the possibility of a vanishing depth while keeping the dispersive 
effects is more difficult, see [15|, |4l| for instance. As for the breaking of waves, 
neglecting the nonhydrostatic and dispersive effects makes the things simpler. 
Indeed, various efficient schemes have been developed to handle the possibility 
of vanishing depth for the NSWE with source terms, relying for instance on 
coordinates transformations [8| , artificial porosity 38| or even variables extrap- 
olations [25|. In a simpler way, it is shown in [27[ that the occurrence of dry 
areas can be naturally handled with a water height positivity preserving finite 
volume scheme, without introducing any numerical trick. Of course, the price 
to pay is the same as above: the dispersive effects are lost. Hence the third 
motivation for this paper: 

tt 3 Propose a simple numerical method that allows at the same time the 
possibility of vanishing depth and dispersive effects. 

The strategy adopted here to handle correctly the three difficulties (t 1-3 
identified above starts from the Green-Naghdi equations. As already said, they 
are very well adapted to tt 1- With a careful choice of the numerical methods, 
they also allow for the possibility of vanishing depth, and thus answer to tt 
3. The main difficulty is thus to handle tt 2 (i.e. wave breaking) with a code 
based on the Green-Naghdi equation. In order to do so, we use a numerical 



scheme that decomposes the hyperbohc and dispersive parts of the equations 
|17| . We also refer to [26[ for a recent numerical analysis of the Green-Naghdi 
equations based on a Godunov type scheme and that provide good results for the 
dam break problem. We use here a second order splitting scheme, we compute 
the approximation U"'~^^ = ((^"+i^ ]/"+i) at time (n + l)St in terms of the 
approximation f7" at time nSt by solving 

U''+^ ^ Si{St/2)S2{St)Si{St/2)U", 

where Si{-) is the solution operator associated to the NSWE (I? = in ([T])) and 
S2(-) the solution operator associated to the dispersive part of the equations 
(keeping only the time derivatives and the r.h.s. of ([T])). For the numerical 
computation of S'i(-), we use a high order, robust and well-balanced finite vol- 
ume method, based on a relaxation approach [3|. This method is known to be 
computationally cheap and very eflficient to handle wave breaking and presents 
another interesting feature for our purposes: it allows the localisation of the 
shocks. In the vicinity of these shocks (or bores to use the physical term) the 
derivation of the dispersive components of the Green-Naghdi equation is mean- 
ingless and these terms, which contain third order derivatives, become very 
singular; moreover, it is known [6|, |9[ that the NSWE correctly describe the 
dynamics of the waves near the breaking point. We therefore "skip" the com- 
putation of S2{-) near the shocks detected during the computation of Si{St/2). 
Elsewhere, 5*2 (•) is computed using a finite diflFerence scheme (note that a careful 
mathematical analysis of 52(-) allows considerable simplifications and numerical 
improvements) . 

In Section[21 we present the physical model studied here, namely, the Green- 
Naghdi equations. After giving the formulation of the equations in non-dimen- 
sionalized form in ij2.11 we show in ^2.21 that it is possible to rewrite them in 
a convenient way that does not require the computation of any third order 
derivative of the unknowns C and V, and exploits the regularizing properties 
of the equations. With classical methods, we then turn to derive a family of 
Green-Naghdi equations with improved frequency dispersion in W2.31 depending 
on a parameter a to be chosen. Another, still non-dimensionlized, reformulation 
of the equations (in terms of {h, hV) rather than {(, V)) is then given in ^2.41 
finally, we give a version with dimensions of these equations in ^2.51 

Section 13] is then devoted to the presentation of the numerical scheme. The 
hyperbolic/dispersive splitting is introduced in WS.W We then turn to describe 
the spacial discretization of the hyperbolic and dispersive parts in H3.2I and 
^3.3\ respectively. The time discretization is described in ^^ where the con- 
sequences of our approach for the dispersive properties of the model are also 
studied carefully. We show in particular that it is possible to derive an exact 
formula for the semi-discrete dispersion relation that approaches the exact one 
at order two. We then use this formula to choose the best parameter a for 
the frequency-improved Green-Naghdi equations; this choice differs from the 
classical one based on the exact dispersion relation. 

Finally, we present in Section 2] several numerical validations of our model. 
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Figure 1: Sketch of the domain 

We first consider the case of solitary waves in §4.11 and use it as a validation 
tool for our numerical scheme. We then evaluate the dispersive properties of 
the model by considering the propagation of a periodic and regular wave over 
a flat bottom in ^4.21 this test illustrates the interest of choosing the frequency 
parameter a in terms of the semi-discrete dispersion relation. In ^4.3[ wc focus 
on the reflection of a solitary wave at a wall, while the ability of the model to 
simulate the nonlinear shoaling of solitary waves over regular sloping beaches 
is investigated in ^4.41 At last, the run-up and run-down of a breaking solitary 
wave is studied in 



2. The physical model 

Throughout this paper, we denote by C(i, X) the elevation of the surface with 
respect to its rest state, and by —ho + b{X) a parametrization of the bottom, 
where ho is a reference depth (see Figure [1]). Here X stands for the horizontal 
variables {X = {x,y) for 2D surface waves, and X = x ioi ID surface waves), 
and t is the time variable; we also denote by z the vertical variable. 
If Uhor denotes the horizontal component of the velocity field in the fiuid domain, 
we then define V as 

1 f^ 

V{t,X) = - / Uhor{t,X,z)dz, 

"- J-ha+b 

where /i := /lo + C ~ ^ is the water depth. We thus have V = {u, v) G M^ for 2D 
surface waves, and V = u gM. for ID surface waves. 

Denoting by a the typical amplitude of the waves, by abott the typical amplitude 
of the bottom variations, and by A the order of the wavelength of the wave, it 
is possible to define dimensionless variables and unknowns as 

and 

C T b r, V 



C=-, b= , V 

a abott 



We also define three dimensionless parameters as 

_ _^ _ ^ „ _ atott ^ 

ho V ho 

here e denotes the nonUnearity parameter, /i is the shallowness parameter while 
(3 accounts for the topography variations. 

2.1. The non-dimensionalized Green-Naghdi equations 

According to [l|, |2J] , the Green-Naghdi equations can be written under the 
following non-dimensionalized form (wc omit the tildes for dimensionless quan- 
tities for the sake of clarity) : 



(3) 



9tC + V • (hV) ^ 0, 

(/ + nT[h, b])dtV + VC + s{V ■ V)y + etiQ[h, b]{V) = 0, 

where we still denote by h the non-dimensionalized water depth, 

h = \ + eC, - Ph, 

and the linear operator T[h, b]- and the quadratic form Q[h^ 6](-) are defined for 
all smooth enough R'^-valued function M^ (d = 1, 2 is the surface dimension) by 

T[h,h]W ^ni[Kb]{\I ■W)+ pn2[h,b]{\Ib-W) (4) 

Q[h, b]{W) ^TZi[h, 6](V • {WV ■ W) - 2(V • W)^) + /37^2[/^, 6]((W • V)^^), (5) 

with, for all smooth enough scalar-valued function w, 

TZi[h,b]w = -l-V{h^w)-l3-wVb, (6) 

'R2[h,b]w = ^V(/i2u;)+/3u;V6. (7) 

Notation 2.1. For the sake of clarity, and one no confusion is possible, we 
often write T, Q, 7?,i andTZ2 instead ofT[h,b], Q[h,b], etc. 



Remark 2.2. For practical applications (see for instance in ^4-5\ ), a classical 
quadratic friction term can be added to the right-hand side of the momentum 
equation. It has the following expression: — fefi^^''^ j^\\V\\V , where f is a non- 
dimensional friction coefficient. 

2.2. First reformulation of the equations 
Let us now remark that 

Q{V)^T{{V-\/)V) + Q,{V), 

where Qi{V) only involves second order derivatives of V (while third order 
derivatives appear in Q{V)); indeed, a close look at (j4)) and ([5]) shows that 

Qi{v) = ni{v -(yv -v -{v-v)v)-2{v -vf) 

+/37^2(V•(T^•V)V6) 

= -27^l(aly • d2V^ + (v • vf) + pn2{v ■ iy ■ v)v&). 



with V^ = (— V2,Fi)^. The fact that this expression does not involve third 
order derivatives is of great interest for the numerical applications. 
We have thus obtained the following equivalent formulation of the Green- Naghdi 
equations ([3]): 

r a^C + V • (hv) = 0, 

\ (/ + i-iT)dtV + e{I + fiT){V ■ W)V + VC + efiQi{V) = 0, ^ ' 

with h = 1 + eC — /3b, and where the quadratic form Qi is given by 

Qi[h, b]{V) = -2ni{diV ■ d2V^ + (V • Vf) + p'R2{V ■ {V ■ V)V6) (9) 
(the linear operators T, TZx and 7?.2 being defined in (|3]), ([H]) and ([7]). 
2.3. Green-Naghdi equations with improved dispersive properties 



It is classical |42|, |28| or [15| that the frequency dispersion of ([5]) can be im- 
proved by adding some terms of order 0(/j,^) to the momentum equation. Since 
this equation is already precise up to terms of order 0(/z^), this manipulation 
does not affect the precision of the model. Such a manipulation is also per- 
formed in (26| but with the goal of working with potential variables rather than 
improving the frequency dispersion. 
The first step consists in noticing that, from the second equation in ([5]), one has 

dtV = -VC - e{V ■ V)V + 0{p), 

and therefore, for any parameter a € M, 

dtV = adtV - (1 - a)(VC + e{V ■ \I)V) + 0(/^). 

Replacing dtV by this expression in ([8]) and dropping the 0{fj,^) terms yields 
the following Green-Naghdi equations with improved frequency dispersion, 

{I + fiaT)dtV + e{I + naT){V ■ V)V (10) 

+ (/ - fi{l - a)r)VC + sfiQiiV) = 0. 

Of course, © corresponds to a particular case of ((TU)) with a = 1. The interest 
of working with (jlO[) is that it allows to improve the dispersive properties of 
the model by minimizing - thanks to the parameter a - the phase velocity error 
(see l2.5p . In [1^, a three-parameter family of formally equivalent Green-Naghdi 
equations is derived yielding further improvements of the dispersive properties. 
For the sake of simplicity, we stick here to the one-parameter family of Green- 
Naghdi systems (fTOj) . 

2.4. Reformulation in terms of the {h, hV) variables 

The Green-Naghdi equations with improved dispersion ()10p are stated as two 
evolution equations for C and V. It is possible to give an equivalent formulation 



as a system of two evolution equations on h and hV, as shown in this section. 
For the first equation, one just has to remark that edtC = dth, so that 

dth + eV -{hV) ==0. 

For the second equation, we first use this identity to remark that hdtV = 
dt{hV) + eV • {hV)V. Multiplying the second equation of pH)) by h, and using 
the identity 

V ■{hV(g)V)^V ■ {hV)V + h{V ■ V)l/, 

we thus get 

(/ + tiahTj^)dt{hV) + e{I + ^iahTj^)\' ■ [hV (g> V) 
+ {I - fi{l - a)hT{)hVC + e^lhQl{V) = 0. 

The Green-Naghdi equations with improved dispersion can therefore be written 
in (h, hV) variables as 

dth + eV -{hV) ==0, 

(/ + fiahrhdtihV) + e{I + fiahrhv ■ {hV ® V) ^^^ 

+ {I- M(f - a)hT-r)hWC + SfihQiiV) = 0. 
h 

The second equation of PT|) requires the computation of third order derivatives 
of ( that can be numerically stiff. It is however possible to show that these 
terms can be factorized by / + fxahT^ , up to a term involving only a first order 
derivative of C,: 

(I - ^l{l - a)hrT)h\/C = -/iVC + ^^^(/ + fiahrhhX^C. 
h a a h 

The equations PT|) can therefore be reformulated as 

dth + eV -{hV) = 0, 

a - 1 
dt{hV)+eV ■{hV(S>V) + hVC, (12) 

1 1 '^ 
+ (/ + ^ahT-r)'^[-hVC + e^JihQl{V)] = 0. 
n a 

This formulation does not require the computation of any third-order derivative, 
allowing for more robust numerical computations, especially when the wave 
becomes steeper. 

2.5. Dimensionalized equations 

Going back to variables with dimension, the system of equations (|12p reads 

f dth + V ■ (hV) = 0, 

dt (hV) + ^^ghVC + W ■ihV(E)V) (13) 

11 
+ (/ + ahr-)-^[-ghVC + hQi{V)] = 0, 
h a 
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Linear phase and group velocity errors 




Local optimal values of a 



O^opt 




Figure 2: Top: linear phase velocity (black) and group velocity (grey) errors for a = 1.159 (full 
line, global optimal value) and o = 1 (dashed line, original model). Bottom: local optimal 
values of a for kh in [0; 4] . 



where the dimcnsionahzed version of the operators T and Qi correspond to (|3]), 
(O, (O and (O with /3 = 1, and where h now stands for the water height with 
dimensions, 

h = ho + C-b. 



Looking at the hnearization of (J13p around the rest state h = ha, V ^ 0, and 
flat bottom 6 = 0, one derives the dispersion relation associated to (IT51) . It is 
found by looking for plane wave solutions of the form (h. fey )e'^^'^~"*^ to the 
linearized equations, and consists of two branches parametrized by i-Oa,±{-)-, 



WQ,±(k) = ±|k 




l + (a-l)(|k|/io)V3 



(14) 



As mentioned before, the parameter a can be helpful in optimizing the dis- 
persive properties of the original model {a = 1), namely the linear phase and 
group velocities denoted by Cqj^ and Cqj^. By adjusting a, wc can minimize the 
error relative to the reference phase and group velocities Cg and C| coming from 
Stokes linear theory. A classical approach consists in minimizing the averaged 
error over some range kho & \0;K] (see [ll|, [29[ and [15| for further details). 
Here, minimizing the weighteqj averaged error over the range kho G [0; 4] yields 



^The squared relative error is weighted by 1/kho to keep the errors to a minimum for low 



the global optimal value a = 1.159, that is adopted, unless stated otherwise, 
throughout this papeo In Figure [2] (top), the ratios Cqj^/C^ and Cqj^/Cq 
are plotted against the relative water depth khg for a = 1 (original model) and 
a ~ 1.159 (optimized model). 

The global optimal value a = 1.159 is especially well-suited when consider- 
ing irregular waves or regular waves over uneven bottoms, i.e. when multiple - 
or not easily predictable - wavelengths are involved. However, when considering 
monochromatic waves over flat bottoms, i.e. when only one wavelength is in- 
volved, a = 1.159 is not optimal anymore. In this particular case, an alternative 
approach consists in minimizing the error on the phase velocity for a specific 
value of kho, where k corresponds to the wavcnumbcr of the considered wave. 
For any discrete value of kho, one easily computes the corresponding optimal 
value of a, denoted by OLoptikho) and refcred to as local optimal value. In Figure 
m (bottom), ttopt is plotted against fc/io, for fc/iQ G [0;4]. 

3. Numerical methods 

We propose here to take advantage of the previous reformulation. First, this 
new formulation is well-suited for a splitting approach separating the hyperbolic 
and the dispersive part of the equations P^ . We present our splitting scheme 
in ij3.1l we then show in Ji3.2l and ij3.3l how we treat respectively the hyperbolic 
and dispersive parts of the equations. 

3.1. The splitting scheme 

We decompose the solution operator S{-) associated to P^ at each time 
step by the second order splitting scheme 

S{5t) = Si{5t/2)S2{5t)Si{5t/2), (15) 

where ^i and S2 arc respectively associated to the hyperbolic and dispersive 
parts of the Grecn-Naghdi equations (fT3|) . More precisely: 

• Si (t) is the solution operator associated to NSWE 

dth + y ■ {hV) = 0, 

1 o (16) 

dtihV)+\'{-gh^) + \' ■{hV(E)V) = -ghVb. 

• S2 {t) is the solution operator associated to the remaining (dispersive) part 
of the equations, 

J dth = 0, 

I dt{hV)--ghWC+{I + ahTh'^[-ghWC + hQi{V)] = 0. 
^ a h a' 

(17) 



wavenumbers. 

•^The dispersive correction used in llSf slightly differs from here, yielding a different defini- 
tion of a: denoting by a the parameter used in |15| . the correspondence is given by 5 = "~ . 



10 



Remark 3.1. From this point, we only consider ID surface waves. The nu- 
merical implementation of our scheme for 2D surface waves is left for future 
work. 

Remark 3.2. The friction term (see Remark \2.S^) . when used, is included in 

S2{t). 

Taking advantage of the hyperbolic structure of the NSWE, Si{t) is com- 
puted using a finite-volume approach, as described in the next subsection. As 
far as the operator 5*2 (t) is concerned, we use a finite-difference approach, as 
shown in 



Such a mixed finite-volume finite-difference method implies to work both 
on cell-averaged and nodal values for each unknown. We use the following 
notations: 

Notation 3.3. - The numerical one- dimensional domain D, is uniformly divided 
into N cells (Ci)i<i<jv such that Ci = [xi-i,Xi\, where {xi)o<i<N are the iV+ 1 
nodes of the regular grid. 

- We denote by Sx the cell size and by St the time step. 

- We write w" the nodal value of w at the i*^ node Xi and at time tn = nSt ■ 

- We denote by wf the averaged value of w on the z*'' cell Ci at time tn = nSf. 

The choice of a finite difference method for solving S2 also entails to switch 
between the cell-averaged values and the nodal values of each unknown, in a 
suitable way that preserves the global spatial order of the scheme. Using classical 
fourth-order Taylor expansions, we easily recover the following relations that 
allow to switch between the finite volume unknowns (w")i<i<Ar and the finite 
difference unknowns (w")o<i<iv at each time step: 

1 2 1 1 

-u;,_i + -Wi + -Wi+i = -{ui + u,+i) H- 0{5i), 0<i<N, (18) 



-^u;,_2 + ^w,_i + —w, - ^w,+i + 0{St), l<t<N, (19) 



and 



with adaptations at the boundaries following the method presented in ^3.51 
We can easily check that ([T8| . (|T9)) preserves the steady state at rest, and that 
these formulae are precise up to 0((5^) terms, thus preserving the global order 
of the scheme. 

3.2. Spacial discretization of the hyperbolic component S'i(-) 

When specified in one space dimension, the system under consideration reads 
as follows: 

dth + dxihu) = 0, , . 

dt{hu) + dx{hu^ + gh^/2) = -ghd^b ^ ^ 



fl 



For the sake of simplicity in the notations, it is convenient to rewrite the system 
(|20|) in the following condensed form: 



9tw + a,f(w)=S(w,&), (21) 

with 

(22) 
where w : R x R+ ^ il is the state vector in conservative variables and f (w) : 
ri — >■ R^ stands for the flux function. The convex set il of the admissible states 
is defined by 

n = {vf e R^; h > 0,u e R} . 

Considering numerical approximations of system ()2ip . wc seek a numerical 
scheixie that provides stable simulations of the processes occurring in surf and 
swash areas, with a precise control of the spurious effects induced by numerical 
dissipation and dispersion. Moreover, the scheme should be able to handle the 
complex interactions between waves and topography, including the preservation 
of motionless steady states: 

u — 0, h + b — cste. 

In this way, we use a low-dissipation and well-balanced extension of the robust 
finite volume scheme introduced in [3|. The main features of the first order 
scheme are recalled in §3.2.11 its higher-order and well-balanced extension pre- 
sented respectively in t^3.2.2l and tj3.2.3l 

3.2.1. First order finite-volume scheme for the homogeneous system 
The homogeneous NSWE associated with (PO)) . given by 

9tw + a,f (w) = 0, (23) 

is known to be hyperbolic over 51. As a consequence, the solutions may develop 
shock discontinuities. In order to rule out the unphysical solutions, the system 
(j23p must be supplemented by entropy inequalities (see for instance [7| and 
references therein). 

The spatial discretization of the homogeneous system ((23|) can be recast under 
the following classical semi-discrete finite-volume formalism: 

— Wi(t) + — ff(wi,W,;+i) - f (w,;_i,Wj)j = 

where f is a numerical flux function based on a conservative flux consistent with 
the homogeneous NSWE. For the numerical validations shown in SHI we use the 
numerical flux issued from the relaxation approach introduced in [3| . 

Remark 3.4. The robustness of this finite volume scheme for the homogeneous 
NSWE is shown in |j/, where the detailed study of the relaxation approach is 
performed. 
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3.2.2. A robust high-order extension 

To reduce both numerical dissipation and dispersion within the hyperbohc 
component S'i(-), high order reconstructed states at each interface have to be 



considered. Following the classical MUSCL approach [37[, we consider the mod- 
ified scheme: 

|w,(t) + 1 (f (w^,,wr+i,) - f (wr_i,.,w^O) - 0' (24) 

where w"; and w"^ are high-order interpolated values of the cell-averaged so- 
lution, respectively at the left and right interfaces of the cell Q. The low 
dissipation reconstruction proposed in [l^l is used. Considering a cell Q, and 
the corresponding constant value /i", we introduce linear reconstructed left and 
right values /i"; and /i"^ as follows: 

-hl^ = -h^ + hK:^^ and 1^, = ^-]^^. (25) 

The corresponding gradients are built following the five points stencil: 

ihi, = {\~^m+, - hi) + HK - hu) 

+ e{-hU + S'K - 3/ir+i + K+2) (26) 

ih,, = {i~vm - hu) + ^(/^r+i - m 

+ ei-h^2 + 3K-i~^ + K,i) (27) 

and the coefficients v, ^^ and ^'^ are set respectively to i, —j^ and —j^, leading 
to better dissipation and dispersion properties in the truncature error. 

When the generation of shock waves occurs during computation, the previ- 
ous interpolation has to be embedded into a limitation procedure to keep the 
scheme non oscillatory and positive. We suggest to use a three-entry limitation, 
especially designed to generate a positive scheme of higher possible order far 
from extrema and discontinuities. Scheme (j24p thus becomes 

dt-^^^'^ + l 

The limited high-order reconstructed values are defined, considering for instance 
the water height /i, as 

"^K.^K + huhn and "^hl.^hl-^Luin- (29) 

To define Li_r(^") and Li,i(h"), we use the following limiter: 

L(u V w) ^ I ^ iiuv <Q, 

^ ' ' ' 1 sign(u) min(2|u|, 2|t;|, w) otherwise. ^ ' 
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w,(t) + - (f (^w^:,, ^wr+i,) - f (^wr_i,,„ ^wi;,) = o. (28) 



Relying on ([50)) . wc then define the limiting process as 

Li^rih''''') = L{Sh^' ,6h"' ,Sh'-^) and Li,z(/i") = L((5/i"' , (5/i"' ,Sh2i), 

where 6h.^ ' = hf^i — hf and Sh^ ' = hf — hf_i are upstream and downstream 
variations, and Sh^ ^ and 8h^ ; taken from ([25]) and (P7|) . 

Such limited high order reconstructions must also be performed for the other 
conservative variable hu. 

Remark 3.5. It is straightforward that when the considered conservative vari- 
able is smooth enough, this limiter preserves the high order accuracy of the 
reconstructions 1126]) and \21l^ . In addition, this high order reconstruction and 
the limitation process can easily be extended to non-uniform meshes. 

Remark 3.6. The robustness ofthe resulting high order relaxation scheme can 
be proved following the lines of |j/. 

3.2.3. Well-balancing for steady states 

We finally introduce a well-balanced discretization of the topography source 
term. Scheme psp is embedded within a hydrostatic reconstruction step [7|. 

To achieve both well-balancing and high order accuracy requirements, we 
have to consider not only high order reconstructions of the conservative vari- 
ables, as done in ij3.2.21 but also of the surface elevation C,. The resulting finite 
volume scheme is able to preserve both motionless steady states and water height 
positivity. The reader is referred to [7| for a detailed study of the hydrostatic 
reconstruction method, including robustness, stability and semi-discrete entropy 
inequality results. 

3.3. Spacial discretization of the dispersive component S2{-) 

The system corresponding to the operator 15*2 (•) writes in one dimension 

J dth = 0, 

1 dt{hu)~-ghd^C+il + ahTh-'[-ghd,C + hQ,{u)] = ^^^^ 

where the operators T and Qi are explicitly given by 

Tw = -^dlw - hd.hd.w + (a,C5x6 + ^dlb)w, (32) 

and 

Qi{u) ^ 2hd^ih+^)id^u)'^+^h^d^udlu+hd^bud^u+{d^(:dlb+^dlb)u^. (33) 

As specified in §3.11 the system ([31]) is solved at each time step using a classi- 
cal finite-difference technique. The spatial derivatives are discretized using the 



14 



following fourth-order formulae: 

iSxW)^ = -—-—{-Wi+2+&Wt+i-SWt-l+Wi-2), 

{d'^w)i = — -— (-Wi+2 + 16wi+i - 30wi + 16wi_i - Wi_2), 

{d^w)i = — -3(-mi+3 + 8it;j+2 - 13mj+i + 13w.,_i - 8wj_2 + w,_3). 

Boundary conditions are imposed using the method presented in Wd.5\ 

3.4. Time discretization and dispersive properties 

3.4.1. Time discretization 

As far as time discretization is concerned, we choose to use explicit methods. 
The systems corresponding to Si and ^2 are integrated in time using a classical 
fourth-order Runge-Kutta scheme. 

3.4.2. Dispersive properties 

We now turn to investigate the dispersive properties of our numerical scheme. 
Since the main originality of this approach is the splitting in time of the hyper- 
bolic and dispersive parts, we consider here the scmi-discrctized in time version 
of our numerical scheme. An extension to the fully discrctized scheme is of 
course possible, but extremely technical, and would not bring any significant 
insight on the dispersive properties of the hyperbolic/dispersive splitting. 

We recall that the dispersion relation associated to the Green-Naghdi (with 
improved frequency dispersion) equations ([13]) is given by (fT4)) or, for \D surface 
waves. 



l + a{kha)^/i 

The dispersion relation corresponding to our semi-discretized (in time) splitting 
scheme is given by the following proposition. 

Proposition 1. The dispersion relation associated to the semi-discretized 
scheme i fJ5)) . jjg)). jj?] ) is given by 

Remark 3.7. The proposition above shows that the semi-discretized dispersion 
relation approaches the exact dispersion relation of the Green-Naghdi equations 
il3\) at order 2 in St. An additional information is that the 0{S^) error made 
by the splitting scheme is always real. Therefore, the numerical errors are of 
dispersive type and there is no linear instability induced by the splitting. 
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Remark 3.8. Since the main error in the dispersive relation is of dispersive 
type, it is natural to try to remove it with techniques inspired by the classical 
Lax- Wendroff scheme. This is possible, but this does not yield better results 
than the 6t-optimized choice of the frequency parameter a (see below), which is 
a much simpler method. 

Proof. For the sake of clarity, we still denote by 5'i(-) and 5*2 (•) the solution 
operators associated to the senii-discretized version of the linearization of (|16p 
and ([T7)) around the rest state (and flat bottoms). 
Step 1. We show here that 

V 9hoa2[Ot)ik ai(6t) J 
with 

MSt) = 1 + ${-9hok') + ^{-ghok^f, a^iSt) = S^ + f (-g/^ofc^). 
z z4 D 



Since the linearization of ((16)) around the rest state and flat bottom can be 
written in compact form as 

dtw + A{d,)w = 0, with Aid.:,) = ( ^^H^^ ^" 

the quantity S'i(i5t)(w°e*'^^) corresponding to the RK4 time discretization is 
given by 

Sii6t)iw°e'''^) = (l + 6tA{zk) + ^Aitkf + f ^(ifc)^ + ^A{ik)*)w''e'''\ 
\ 2 24 / 

A simple computation thus yields the result. 
Step 2. We show here that 

^ Tn>2 c /X \r,,,0„ikx\ _ ( 1 \ ^^Si ikx 



Vw>^ e M^ 52(^.)(w^e-^) = ^ ^,^^^^,^ I ) w^e 

with 

[kh^Y/i 
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l + a(fc/io)V3' 

Since the linearization of (J17p around the rest state and flat bottom can be 
written in compact form as 

9,W+i3(9.)w = 0, with B{d.^) = (^ _^^^^^ _ .^2g2)-l(_ 1/^252)9^ 

the quantity 52((5t)(w°e*'^^) corresponding to the RK4 time discretization is 
given by 

^2(^t)(w"e*''") = (1 + <5tB(^fc))w°e'*^^ 

16 



where we used the fact that B{ik)'^ = 0. The resuh follows directly. 
Step 3. By a direct computation, we get that 

Vw" G R^, Si{6t/2)S2{St)Sii6t/2){w°e'''-) = {I + 5tM)w°e'''\ 

with M = {mij)i<ij<2 given by 

b 2 

™2i = igho{l+^)k~ 1^-^(1 + l^)kH^ + 0{6t). 

2 

Step 4. End of the proof. Wc deduce from the previous steps that if -s^o^tkx-ujt 
is a plane wave solution for the semi-discretized scheme, then one has 

e-*"'5*w° = (/ + (5tAf)w°, 

and - — J— ^ is therefore an eigenvalue of M. After some simple computations, 
we thus get 

— iuSt _ 1 

—, -A±, (34) 



-zuj^Mk) - ^^^.+ (fc)'5t + ^c.„,±(fe)^(4 - JY^)St + OiS^). (35) 



with 



By identifying the Taylor expansion of the left-hand-side of ([34]) with (|35|) , we 
deduce that 

^ = w„.± + ^^"^±(^)'^r^:^'^' + '^(^t), (36) 

and the result follows. D 

Starting from the previous expression and dropping the 0{Sf) term, one 
easily obtains the semi-discrete linear phase and group velocities CQj^{5t) and 
CQj^{6t), and computes the semi-discrete error relative to the reference velocities 
Cg and Cg. Obviously, the global value a = 1.159 is no longer optimal with 
the additional 0{S^) term, and we need to compute new optimal values of a 
that depend on St . As in ^2.51 and for discrete values of the non-dimensional 
time step 5t := VT^^tj "^6 look for 1) the global optimal value aopt {St ) over 
the range [0;3], and 2) the local optimal values aoptiSt^kho) for some discrete 
values of kho- 

Results are gathered in Figure [3j the top figure plots the global optimal 
value aopt against St , while the bottom figure plots the local optimal values aopt 
against St, each curve corresponding to a discrete value of kho. Wc point out 
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Global optimal values of a 
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Figure 3: Top: global optimal values of a against St- Bottom: local optimal values of a 
against St for kho = 7r/4 (full grey line), khg = ir/2 (dotted line), khg = 3it/4 (dash-dotted 
line), kho = "" (dashed line), kho = 57r/4 (full black line). 



that in both approaches, the computed optimal value of a was sometimes found 
lower than 1, for instance when St G [0.26; 0.42] for the global optimal value aopt- 
Since taking a < 1 induces some high-frequency instabilities, the optimal value 
of a has been taken equaljbo 1 in such cases. However, it is worth remarking 
that for these problematic (5t-regions, the model that provides the best dispersive 
properties - among the stable ones - is the original Green-Naghdi model. 

We finally refer to ^4.21 for numerical simulations showing the consequences 
of our choice to optimize a taking into account the dispersive effects of the time 
discretization. 



3.5. Boundary conditions 

The boun dary conditions for the hyperbolic part Si of the splitting are 
treated as in |27| . More precisely, as the simulations shown in this work do 
not require complex Riemann invariants based inflow, outflow or absorbing con- 
ditions, we simply introduce "ghosts cells" respectively at left and right bound- 
aries of the domain, and suitable relations are imposed on the cell-averaged 
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quantities : 

• w_fc+i = WN^k+i and WN+k ~ Wk, k > 1, for periodic conditions on the 
left and right boundaries, 

• W-k+i = Wk and iJUN+k = WN-k+i, k > 1, for homogeneous Neumann 
conditions on the left and right boundaries, 

• w^k+i — ~Wk and WN+k ~ ^wpf^k+i, k > 1, for homogeneous Dirichlet 
conditions on the left and right boundaries. 

For the dispersive part S2 of the splitting, the boundary conditions are simply 
imposed by reflecting - periodically for periodic conditions, evenly for Neumann 
conditions and oddly for Dirichlet conditions - the coefficients associated to 
stencil points that are located outside of the domain. The advantage of this 
method is to avoid introducing decentered formulae at the boundaries, while 
maintaining a regular structure in the discretized model: 

• W-k = WN-k and WN+k-i ~ Wk-i, k > 1, for periodic conditions on the 
left and right boundaries, 

• w^k ~ Wk and lOiv+fe = w^-k, k > 1, for homogeneous Neumann condi- 
tions on the left and right boundaries, 

• W-k ~ —Wk and WN+k ~ —wjy-k, k > 1, for homogeneous Dirichlet 
conditions on the left and right boundaries. 

Solid wall effects on the left or right boundary can be easily reproduced by 
imposing an homogeneous Neumann condition on h and h, and an homogeneous 
Dirichlet condition on hu and hu, with the previous methods. 
When the water depth vanishes, a small routine is applied to ensure stability on 
the results: on each cell, if h is smaller than some threshold e then we impose 
the values h = e and v = 0. 

3.6. Wave breaking 

In order to handle wave breaking, we switch from the Green-Naghdi equa- 
tions to the NSWE, locally in time and space, by skipping the dispersive step 
5*2 ((5t) when the wave is ready to break. In this way, we only solve the hyperbolic 
part of the equations for the wave fronts, and the breaking wave dissipation is 
represented by shock energy dissipation (see also [6[, [27[ and [9[). 
To determine where to suppress the dispersive step at each time step, we use 
the first half-time step Si of the time-splitting as a predictor to assess the local 
energy dissipation, given by 

D, = -{dt£ + d^T), (37) 

with £ = ^{hu^ + g(^) the energy density and T— phu{^ + gC) the energy 
flux density. This dissipation is close to zero in regular wave regions, and forms 
a peak when shocks are appearing. We can then easily locate the eventual 
breaking wave fronts at each time step, and skip the dispersive step only at the 
wave fronts. 
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4. Numerical validation 

4.I. Propagation of a solitary wave 

It is known that for horizontal bottoms, the Green-Naghdi model with a = 1 
have an exact solitary wave solution given by 

h(x,t) = ho + a scch (^(a; — ct)), 

"^"'^^ = <^-/^)' (38) 

c=^/g{h + H), 



2hoVho + H' 



This family of solutions can be used as a validation tool for our present numer- 
ical scheme. We successively consider the propagation of two solitary waves of 
different relative amplitude a/h^, on a 30 ?Ti long domain with a constant depth 
ho = 0.5 m. The considered relative amplitudes are a/ho ~ 0.05 for the first 
solitary wave, and a/ho ~ 0.2 for the second one. Periodic conditions are im- 
posed on each boundary, and the initial surface and velocity profiles are centered 
at the middle of the domain. 

In order to assess the convergence of our numerical scheme, the numerical so- 
lution is computed for several time steps St and cell sizes Sx-, over a sufficient 
duration T — 3s. Starting with Sx = ^rn and St = Sx/^/gho = 0.45s, we 
successively divide the time step by two, while keeping the CFL equal to 1. For 
each computation and each discrete time tn = nSt, the relative errors Ep and 
-E" on the free surface elevation and the averaged velocity are computed using 
the discrete L°° norm ||.||oo : 

j-^ri W'^nuni '^solWoo j^n W^nuni ^so/||oo 

i^r = —rr, ; — r. ; A, = r. r. 



\\hsol — ho\ 



\Usol\ 



where {hnum,u„um) are the numerical solutions and (CsohUsoi) denotes the an- 
alytical ones coming from (|38|) . 

Results are gathered in Figure 01 where maxEQ is plotted against St, for the 
two considered relative amplitudes a/ho = 0-05 and a/ho ~ 0.2. In both cases, 
the convergence of our numerical scheme is clearly demonstrated. Furthermore, 
computing a linear regression on all points yields a slope equal to 1.91 for the 
first case and 1.83 for the second one. This result is coherent since the global 
order of our scheme is obviously limited by the order of the splitting method 
used here, which is of order two. 

4-2. Propagation of periodic and regular nonlinear waves 

In this test, we want to evaluate the dispersive properties of the model, along 
with the optimisation on the semi-discrete dispersion relation proposed in ^13.4.21 
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Figure 4: Propagation of a solitary wave over a flat bottom, maximum of the relative error 
on the free surface elevation. Case a/ho = 0.05 in black, case a/ho = 0.2 in grey. 



We consider the propagation of two-dimensional periodic and regular nonlinear 
waves, without change of form, over a flat bottom. For this situation, numerical 
reference solutions can be obtained by the so-called stream function method (see 
[33[). Unlike analytical wave theories (such as Stokes or cnoidal wave theories), 
this numerical approach is applicable whatever the shallowness and nonlinearity 
parameters /x and e are, and very accurate solutions can be obtained with a 
high number of terms in the Fourier series. These reference solutions are here 
obtained with the software Stream_HT, implemented by Benoit et al. [5[. 

We consider a domain which covers one wave-length {L = X = 2 m), and a still 
water depth ho = 1 to, so that the relative water depth is khg = tt w 3.14. The 
wave amplitude is a = 0.01 ?7i, so that the nonlinearity parameter is a//io = 0.01. 
These conditions correspond to very dispersive and weakly nonlinear waves. 

The domain is discretized with 50 cells {Sx — 0.04 m), and a time-step 5t ~ 0.03 s 
(corresponding to a Courant number Cr = 2.3) is used during the simulations. 
Periodic conditions are imposed at the two lateral boundaries. The period 
computed with the stream function approach (at order 20) is T = 1.133 s and the 
solutions obtained for the the water height and the averaged velocity are imposed 
as initial conditions. Numerical integration is performed over a duration of 25 T. 

Two different values of a are considered: a ~ 1.16, corresponding to the local 
optimal value computed for kho = tt as in ^2.5[ and a ~ 1.153, correspond- 
ing to the local optimal value computed for kho — n and 5t ~ 0.094s, as in 
§3.4.21 The local minimization approach is prefcred since the considered wave 
is monochromatic. 

Results after 25 periods are gathered in Figure [5l where the water height 
(left) and the averaged velocity (right) are plotted and compared to the reference 
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Figure 5: Numerical results at t = 25T, model with a = 1.16 (top) and a = 1.153 (bottom) 
in dashed line, reference solution in full line. 





Relative error on the 
wave amplitude 


Relative error on the 
wave celerity 


a 


= 1.16 


1.8 10-2 


5.10-3 


a - 


= 1.153 


1.7.10-2 


8.10--^ 



Table 1: Relative errors on the wave amplitude and the wave celerity between the model 
results and the reference solution (see Figure [5| at t = 25 T. 



solutions (propagating at constant speed and without change of form). The 
results obtained with a = 1.153 are seen to be in excellent agreement with the 
reference solution, whereas the ones obtained with a = 1.16 are less satisfying. 
However, we point out that these latter provide an overall good agreement with 
the reference solution, as shown in Table 1 where the relative error on the wave 
amplitude at t = 25 T and the relative error on the wave celerity (using the 
phase shift between the solutions at t = 25 T) have been computed. To sum up, 
this test clearly demonstrates 1) the ability of our model to handle intermediate 
or even deep water waves (even for non-optimal values of a), and 2) the interest 
of using the optimisation on the semi-discrete dispersion relation proposed in 
ij3.4.2l when dispersive waves are involved. 

J^.3. Reflection of a solitary wave at a wall 

In this test, we compare numerical results with experimental data taken 
from [39[ , for the propagation and reflexion of a solitary wave against a vertical 
wall. The depth profile, a sloping beach of 1:50 terminated by the wall located 
at X = 771, is depicted in Figure \6\ The aim of this test is to study the full 
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Figure 6: Topography layout for the reflection of a sohtary wave at a wall. 



H =0.07, 500 paints, oi 




10 1? 14 16 11 




Figure 7: Reflection of a solitary wave against a vertical wall. Time series of the free surface 
at a; = 17.75 m for a = 0.07 (left) and a = 0.12 (right): comparison between numerical results 
(blue line) and experimental data (green line). 



reflection of a non-breaking solitary wave propagating above a regular sloping 
beach, before reaching a vertical solid wall. 

The spatial domain is 60 m long, the initial solitary wave is centred at a; = 50 m 
and is propagating from right to left. The still water depth is ho = 0.7 m. The 
boundary condition at right is open, as there is no inflow. At the left boundary 
we use solid walls fully reflective conditions. 

Two runs are performed with two different initial solitary wave amplitudes, given 
in terms of relative amplitude a/ho = 0.1 and a/ho = 0.174. The computational 
domain is discretized using 500 cells and a time step St = 0.05 s is used. 

Numerical results are shown as time series of the surface elevation, at a 
location near the solid wall {x = 17.75 tti). Experimental data are compared 
with numerical results on Figure [T] We can observe the two expected peaks 
corresponding respectively to the incident and reflected waves. We can observe 
a very accurate matching between simulation and experimental data, even for 
the second simulation which involves a more complex non-linear propagation. 
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Incident wave amplitude: 


ao//io = 


0.096 






Gauge location (m) 
(relative to the shoreline) 
Relative amplitude error (%) 


2.430 
-1.6 


2.215 
-2.5 


1.960 
-5.5 


1.740 
-7.1 


1.502 
-10.9 


Incident wave 


amplitude: 


ao/^o = 


0.298 






Gauge location (m) 
(relative to the shoreline) 
Relative amplitude error (%) 


3.980 
1.2 


3.765 
-0.5 


3.510 
0.2 


3.290 
-0.2 


3.052 
0.04 


Incident wave 


amplitude: 


ao/^o = 


0.456 






Gauge location (m) 
(relative to the shoreline) 
Relative amplitude error (%) 


4.910 
3.6 


4.695 
-0.3 


4.440 
1.1 


4.220 
0.5 


3.982 
2.2 


Incident wave 


amplitude: 


ao//io = 


0.534 






Gauge location (m) 
(relative to the shoreline) 
Relative amplitude error (%) 


5.180 
0.03 


4.965 
-0.1 


4.710 
-1.4 


4.490 
-1.7 


4.252 
0.7 



Table 2: Location of wave gauges for solitary waves shoaling on a 1:30 sloped beach, and 
relative error between the computed and measured wave amplitudes at each gauge. 

4.4- Nonlinear shoaling of solitary waves propagating over a beach 

We investigate in this test the ability of the scheme to simulate the nonlinear 
shoaling of solitary waves over regular sloping beaches, which is a paramount 
in the study of nearshore propagating waves. This test is based on the ex- 
periments performed at the LEGI, in Grenoble (France) and reported in |20| . 
Solitary waves are generated in a 36 m long wave-flume, following the procedure 
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Free surface displacements were measured at various locations, using wave gau- 
ges located just before breaking. Four solitary waves of different heights are 
generated (see Table [H, in order to account for various nonlinearity effects 
during propagation towards the shore. 

All simulations are performed using 6x = 0.025 m and 6t = 0.016 s. The initial 
water depth is hg = 0.25 m in the horizontal part of the channel. 

Results are shown for each configuration in terms of time-series at the wave 
gauges locations in Figure[51 The relative error between computed and measured 
wave amplitudes is presented in Table [5] The global agreement is good, both 
for the amplitude and shape of the solitary waves. Significant errors can be 
observed for the less nonlinear case {ao/ho = 0.096), but the discrepancies can 
be partly explained by experimental problems, since it can be observed that the 
water surface is not totally at rest before the propagation of the solitary wave . 

4-5. Run-up and run-down of a breaking solitary wave over a planar beach 

This test is based on experiments carried out by Synolakis |36| for an incident 
solitary wave of relative amplitude ao/ho = 0.28, which propagates and breaks 
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Figure 8: Nonlinear shoaling of solitary waves propagating over a beach - Time series of the 
free surface elevation for the solitary wave propagating over the 1:30 sloping beach. ( — ) 
experimental data, ( ) numerical results, with t* = t(g//io)^'^. 



over a planar beach with a slope of 1:19.85. Free surface elevations at different 
times are available thanks to video measurements. 

The still water level in the horizontal part of the beach is ft-o = 0.3 771. The 
simulations are performed using the cell size E^ = 0.08 777 and (5t = 0.02 s. Fric- 
tion effects are important when the water becomes very shallow, as for instance 
in the run-up and run-down stage. To take into account this phenomenon, we 
introduce a quadratic friction term with the friction coefficient / = 0.002 (see 
Remark !^ . 

The comparison between measured and computed waves is presented in Fi- 
gure [HI It shows a good agreement between model predictions and laboratory 
data and illustrates the ability of our model to reproduce shoaling, breaking, 
run-up and run-down, as well as the formation and breaking of the backwash 
bore, and this without any additional treatment. 

5. Conclusion 

In this work, the fully nonlinear and weakly dispersive Green-Naghdi model 
for shallow water waves is considered. The original model is reformulated un- 



25 



0.4 






/ 


0.4 






/ 


0.2 


/ 




/ 


0.2 


; 


'1 


/ 





/ \ 




' 





HHH^ 




' 




0.2 




/ 




-0.2 




/ 




20 


30 


40 


20 


30 


40 






f=21 








f=22 






Figure 9: Comparisons of model predictions ( — ) and experimental snapshots (+) for a break- 
ing solitary wave with non-dimensional initial incident amplitude ao/ho = 0.28, on a 1 : 19.85 
constant slope beach investigated by Synolakis (1987), with t* = t(g/ho)^'^. 



der a more convenient form, and variants with improved frequency dispersion 
depending on a parameter a are derived. 

A hybrid finite-volume and finite-difference method is then implemented, 
embedded in a splitting approach especially designed to describe the wide range 
of phenomena encountered in coastal oceanography. The first component of the 
model, regarded as a set of hyperbolic conservation laws, is discretized using 
an efficient and robust Godunov-like high-order accuracy finite- volume scheme, 
while high-order finite-differences are used for the dispersive part of the model. 

The dispersive properties of the splitted semi-discrete scheme are carefully 
studied and are shown to approach the dispersion relation of the Green- Naghdi 
model at order two. Moreover, we use the explicit formula for the scmi-discre- 
tized dispersion relation to choose the best coefficient a for the frequency- 
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improved GN model. This optimal value depends on the time step St, and 
is shown to provide better results than the standard choice based on the disper- 
sion relation of the mathematical model. This is because this choice takes into 
account the dispersive effects due to the time discretization. 

In a last part, this new scheme is widely validated. Analytical solutions 
for propagating solitary waves are first considered and allow us to study the 
accuracy and convergence properties of this approach. In the following cases, 
numerical results are compared in an extensive way with both experimental 
data and reference solutions. In particular, the propagation and shoaling of 
highly nonlinear waves are successfully described, together with wave break- 
ing and subsequent run-up and back-wash over a slopping beach. This clearly 
demonstrates the validity of this shock-capturing finite- volume based approach 
for dispersive waves, which appears therefore as a promising tool for the study 
of shallow water waves in coastal areas. 

Following the steps of this study, next steps may concerns the derivation of 
dispersion optimized models |12| , the design of numerical sensor for the accurate 
detection of breaking waves and of course two-dimensional extension of the 
numerical scheme to study more realistic cases. 
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